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rS 1 INTRODUCTION 



ABSTRACT 

We examine predictions for the quasar luminosity functions (QLF) and quasar clus- 
tering at high redshift {z ^ 4.75) using MassiveBlack, our new hydrodynamic cosmo- 
logical simulation which includes a self-consistent model for black hole growth and 
feedback. We show that the model reproduces the Sloan QLF within observational 
constraints at z ^ 5. We find that the high-z QLF is consistent with a redshift- 
independent occupation distribution of BHs among dark matter halos (which we pro- 
vide) such that the evolution of the QLF follows that of the halo mass function. The 
sole exception is the bright-end at z = 6 and 7, where BHs in high-mass halos tend 
to be unusually bright due to extended periods of Eddington growth caused by high 
density cold flows into the halo center. We further use these luminosity functions to 
make predictions for the number density of quasars in upcoming surveys, predicting 
there should be ~ 119 ± 28 (~ 87 ± 28) quasars detectable in the F125W band of 
the WIDE (DEEP) fields of the Cosmic Assembly Near-infrared Deep Extragalactic 
Legacy Survey (CANDELS) from z = 5 - 6, ~ 19 ± 7 (~ 18 ± 9) from z = 6-7, 
and ~ 1.7±1.5 (~ 1.5±1.5) from z — 7 — 8. We also investigate quasar clustering, 
finding that the correlation length is fully consistent with current constraints for Sloan 
quasars (ro ~ 17 h~^ Mpc at z = 4 for quasars above rrii — 20.2), and grows slowly 
with redshift up to z = 6 {tq ^ 22 h^^ Mpc). Finally, we note that the quasar cluster- 
ing strength depends weakly on luminosity for low Lbh, but gets stronger at higher 
Lbh as the BHs are found in higher mass halos. 

Key words: quasars: general — galaxies: active — black hole physics — methods: 
numerical — galaxies: haloes 



Malbon e t all l2007l: ICiotti fc Ostrikeij l2007l : ISiiacki et all 
2007; Hop kins et atlbOOTl ). 



Supemiassive black holes are bel ieved to be present at 
the center of most galaxies l|Korniendv fc Richstond 
[l995|) and are found to be correl ated with the 
prop er ties of their ho st ga la xies (IMagorrian et al.l 

[Cebhardt et al.l l2000l : 



19981: iFerrarese fc MerrittI |2000|: 



Tremaine et al.l |2002| : iGraham fc Driveij 



20071 ). These 

correlations provide strong evidence for a link between 
the growth of black holes and the evolution of their host 
galaxies, g enerally attribu t ed to some form of quasar 
feedback (iBurkert fc Silkl l200l|: ICranato et al.1 120041; 



Sazonov et al.ll2004l: ISpringel et af 
20051; iKawata fc GibsonI l2005l: 



2005bl: IChurazov et all 
Di Matteo et al.1 12005" 



Perhaps the most fundamental statistical quan- 
tity in the study of quasars is the number density, 
often characterized as a Quasar Luminosity Function 
QLF). The QLF h as bee n studied observational!: 



La Franca et a l. 2002; Fiore et al 
Barger et a lj [2 003h: .Groom et al 



2005 



Girasuolo et al.l 



20091 : lYencho et all 



20051: 



l2003l; lUeda et al.|[200: 
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l2004l: iLa Franca et all 



Richards et al.l 



2006 



Brown et al.l l2006l: I Silverman et al.1 l2008l : lEbrero et all 

through simula- 



tions, with black 



20091 ). as well as 
holes modeled both 



Bower et al.ll2006l : iBegelman et aLll2006l : iGroton et al.ll2006l : I2OO8I : IBonoU et al.ll2009l ) and directly incorporated in the 



semi-ana ly tically 



JKauffmann fc Hae hneitl I2OO0I: IVolonteri et al.l 12003: 



Wyithc fc Locb 2 00^: iGranato et all |2004| : iMaruUi et all 
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simu l ation jHopkins et al.ll2005bl , l2006bl , l2007l : iMarulli et all 



I2OO9I : iDeGraf et al.ll201ol T 



Table 1. Numerical Parameters 



Another fundamental quantity for quasar populations 
is the strength of quasar clustering, and its evolution with 
redshift. Numerous ob s ervational studies have been done 
(*La Franca et al.1 [l998l; Porciani et alj 2004 Groom et al 



200.5; Shen et al 



20081 : IShen et al. 



2007] : iMvers et al.1 l2007l : Ida Angela et al 



20091: iRoss et al. 



20091 ) 



generally finding 
evidence for increasing clustering amplitude with redshift, 
in agreement with findings from simu lations (|Bonoli et al.l 
120091 : lGrotonll2009l : iDeGraf et al.ll2010l ). Glustering is espe- 
cially significant because one can use it to estimate the typ- 
ical dark matter halos in which quasars are found simply 
by comparing the clustering strength of quasars to that of 
dark matter halos. In particular, the luminosity dependence 
of the clustering (if any) can help determine if the bright 
and faint quasars populate the same halos, suggesting they 
may be the same pop ulation of quasars at differ ent pha ses 
of their lives (see, e.g. iHopkins et aI]l2005bl ldBfl l2006al ) or 
different halo masses, suggesting th ey are fundament ally dif- 
ferent quasar populations (see, e.g. iLidz et al.ll2006r ). 

Thus both QLF and quasar clustering can be used to 
investigate the populations of quasars being observed, how 
they populate dark matter halos, and how the typical lu- 
minosities correlate with hosts. However, investigations at 
high redshift have been extremely difficult. Simulations have 
proved difficult due to the volumes needed to produce high 
redshift quasars in sufficient numbers, limiting the statis- 
tical investigations. Similarly, the difficulty observing such 
distant objects limits the current number of observed high 
redshift quasars such that the z ~ 6 faint-end QLF is de- 
termined by only a few objects (|Willott et al.ll2010br ), and 
high-redshift quasar clustering is limited to lower redshift 
(e.g. 2 ~ 4 bv lShen et al.ll2009l ). However, upcoming surveys 
such as the James Webb Space Telescope (JWST) and the 
Gosmic Assembly Near-infrared Deep Extragalactic Legacy 
Survey (GANDELS) have the potential to drastically im- 
prove the high-z observations, making now an ideal time to 
make predictions as to the statistics of high redshift quasars. 

In this paper we use a new large-scale hydrodynamic 
cosmological simulation to study the earliest supermassive 
black holes, focusing on their luminosity and clustering prop- 
erties. We take advantage of this simulation's large volume 
(it is the largest cosmological simulation which directly mod- 
els the growth and evolution of black holes) to investigate 
the statistical properties of the earliest supermassive black 
holes, focusing on the luminosity function and correlation 
length. We use these quantities to compare with current ob- 
servational measurements, to make predictions for upcom- 
ing surveys, and to investigate the implications of potential 
features which may be detected. 

In Section [2] we describe the numerical sumulation used 
for our analysis. In Section[3]we investigate the Quasar Lu- 
minosity Function (Section 13. l|l and quasar clustering (Sec- 
tion |3]2]) and compare with observations, and our results are 
summarized in Section |4l 



2 METHOD 

In this paper we use a new cosmological hydrodynamic 
simulation of a 533 h~^ Mpc box specifically intended 
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for high-redshift investigations (see Table [TJ. The sim- 
ulation uses the massively parallel cosmolocial TreePM- 
SPH c ode GADGET- 3 (an updated version of GADGET- 
2, see ISpringell 12005^ inc orporating a multi-phase ISM 
model with star formation (|Springel fc Hcrnquist 2003) and 
black hole accretion and feedback ([Springcl ct al. 2005a|; 
iDi Matteo et al.ll2005^ . 

Within our simulation, black holes are modeled as col- 
lisionloss sink particles which form in newly emerging and 
resolved dark matter halos. These halos are found by call- 
ing a friends of friends group finder at regular intervals (in 
time intervals spaced by A log a = log 1.25). Any group 
above a threshold mass of 5 x 10^°h~^ Mq not already 
containing a black hole is provided one by converting its 
densest particle to a sink particle with a seed mass of 
AfaH.socd = 5 X 10^h~^MQ. This seeding prescription is cho- 
sen to reasonably match the expected formation of super- 
massive black holes by gas directly collapsing to BHs with 
Mbh ~ Maeed (e.g. iBromm fc Loebl l2003l : iBegelman et all 
I2OO6I ) or by PopIII stars coUaps i ng to ~ W^Mp BH s at 
z - 30 (|Bromm fc LarsonI |2004| : lYoshida et aLJ l2006i ) fol- 
lowed by sufficient exponential growth to reach Msood by 
the time the host halo reaches ~ 10^" M©. Following in- 
sertion, BHs grow in mass by accretion of surrounding gas 
and by merging with other black holes. Gas is accreted ac- 
cording to Mbh = a '^,l^^"fi)2 (|Hovle fc LvttletonI Il939l : 
iBondi fc Hovlelligil : iBon di 1952i), where p is the local gas 
density, Cs is the local sound speed, v is the velocity of the 
BH relative to the surrounding gas, and a is introduced 
to correct for the reduction of the gas density close to the 
BH due to our effective sub-resolution model for the ISM. 
To allow for the initial rapid BH growth necessary to pro- 
duce sufficiently massive BHs at early time (~ 10^ Mq by 
2: ~ 6) we al low for mildly super-Eddington accretion (con - 
sistent with IVolonteri fc ReesI 120061 : iBegelman et al.l I2OO6I ) , 
but limit it to a maximum of 3 x Msdd to prevent artificially 
high values. 

The BH is assumed to radiate with a bolometric lu- 
minosity proportional to the accretion rate, L — tjMbhc^ 
(Shakura fc Sunvacv 1973), where the radiative efficiency rj 
is fixed to 0.1 throughout the simulation and our analysis. 
To model the expected coupling between the liberated radi- 
ation and the surrounding gas, 5 per cent of the luminosity is 
isotropically deposited to the local black hole kernel as ther- 
mal energy. The 5 per cent value for the coupling factor is 
based on galaxy merger simulations such th at the normaliza- 
tion o f the Mbh — cr relation is reproduced JDi Matteo et al.l 
l2005h . 

The second mode of black hole growth is through merg- 
ers which occur when dark matter halos merge into a single 
halo, such that their black holes fall toward the center of 
the new halo, eventually merging with one another. In cos- 
mological volumes, it is not possible to directly model the 
physics of the infalling BHs at the smallest scales, so a sub- 
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resolution model is used. Since the mergers typically occur 
at the center of a galaxy (i.e. a gas-rich en vironment), we as- 
sume the final coalescence will be rapid (IMakino fc Funatd 
I2OO4I : lEscala et all l2004l : iMaver et al] |2007| ). so we merge 
the BIfs once they are within the spatial resolution of the 
simulation. However, to prevent merging of BHs which are 
rapidly passing one another, mergers are prevented if the 
BHs' velocity relative to one another is too high (compara- 
ble to the local sound speed). 

The model used for black hole creation, accretion and 



feedback has been investiga te d and dis cussed in[S iiacki et al 



2007 ) 



(2008) 



iDi Matteo et al.l 
ICroft et ai] 



200i) : 
20091): 



( 2010l ): iDegraf et all l|201ll ) 



iColberg fc di Mattec 
ISiiacki et al.l \2oM ): 



Mb 



finding it 
relation, the 



iDeGraf et al 

does a good job reproducing the 

total blac k hole mass densit y (jDi Matteo et al.l 120081 ) . 
the QLF (iDeGraf et al.l I2OI0I). and the e xpected black 
hole clustering behavior (jDegraf et al.l I2OIII ). This simple 
model thus appears to model the growth, activity, and 
evolution of supermassive black holes in a cosmological 
context surprisingly well (though the detailed treatement 
of the accretion physics is infe asible for cosmological scale 
simulations). We also n ote that iBooth fc Schavd ( 2003 ) and 
[Johansson et al.l l|2008l ) have adopted a very similar model, 
and have independently investigated the parameter space 
of the reference model of lDi Matteo et al.l (|2008h . as well as 
varying some of the underlying prescriptions. For further 
details on the simulation metho ds and convergence stu dies 
done for similar simulations, see lDi Matteo et al.l (|2008l ). 

Because the simulation saves the complete set of black 
hole properties (mass, accretion rate, position, local gas den- 
sity, sound speed, velocity, and BH velocity relative to local 
gas) for each BH at every timestep, the black hole output 
for such a large simulation is prohibatively dif ficult to ana- 
lyze u sing previous techniques. For this reason, iLopez et al.l 
l|201 J) developed a relational database management system 
specifically for this simulation. A similar strategy has also 
been followed in the analysis of th e Millenium simulation 
l|Lemson Sz Virgo Consortium! 120061 ). In addition to provid- 
ing a substantially more efficient query system for extract- 
ing information, this database is significantly more flexible 
than traditional approaches. For a complete su mmary of the 
datab ase format and its efficiency, please see iLopez et al.l 
(|201ll ). 



3 RESULTS 

3.1 Luminosity Function 

In the left panel of Figure [T] we show the bolometric Quasar 
Luminosity Function for Lbh > IO^^Lq at z = 11 to 2: = 5 
(solid li nes). We also show the observational data com- 
piled bylWfllott et all (|2010ah from SDSS main (diamonds, 
Fan et allbOOfil ). SDSS deep stripe (triangles, Ijiang et all 
20091). and the Caud a- France High-z Quasar Survey (circles, 
Willott et al.ll2010bl), and t hat compiled by Hopkins et al. 
ll2007h (asterisk s; iFanet al.li2001b,a; Bargcr c t al. 2003b;a,; 
Fan et all I200I 12004 ICristiani et al.1 |2004|: iBarger et all 



20051 : iRichards et al.ll2005l : [Silverman et al.ll2005h . We find 
that our simulation is generally consistent with observations, 
though we tend to predict a steeper slope than observations. 



Thus this large simulation volume shows that our model 
is capable of producing Sloan-type BHs at high redshift 
{z = 5,6), and in the co rrect abunda nces (and in fact of 
sufficiently high mass, see lDi Matteo et al. 2011 ). 

To help better understand the quasar populations which 
we are simulating and the dark matter halos in which they 
are found, we have performed a simple fit to characterize 
the relation between BH luminosity and host halo mass. 
To do this, we use a simple Halo Occupation Distribution 
(HOD) model. Because we are only interested in the high- 
luminosity sources (and it is exceedingly rare to find multiple 
high-luminosity objects within a single halo, especially at 
such high redshift), we model the probability that a halo 
of a given mass (Mhos*) will host a BH above a specified 
luminosity (I/BH,cut) as a cumulative log-normal distribution 
(wh ich is often used in galaxy HOD models for (Nccn), see, 
e.g. IZheng et al.ll2005l ) 



/(LBH,c.t|AfHoeO = iU + erfr^(^«-* 



m(Le 



■y/2cr2(LBH,cut) 

where /(I/BH,cut|AfHost) is the fraction of halos with mass 
AfHost which contain a BH above Lbh.cui, and erf is the er- 
ror function I erf (x) = -^ P e~^ dy J [see Table[2]for fitted 

parameters] . We combine this occupation fractio n (using pa- 
ramet ers at 2 = 5) with the halo mass function of lReed et al.l 
((20071) to predict the QLF, shown in Figure [l] (blue dashed 
line). This predicted QLF matches the simulation well at 
z = 5, confirming an accurate fit. We note that our HOD ap- 
proach may slightly underpredict the low luminosities since 
the HOD parameters are fit for bright (> 10^ L©) quasars, 
and because at low luminosities, satellite BHs (which are 
not included in our simple HOD) may become significant. 
For bright objects, however, this approach is very accurate. 

We also plot the QLF for higher redshifts using the 
z = 5 HOD fit (dashed lines), and find that it does a surpris- 
ingly good job of predicting the higher redshift QLF, with 
the expection of the luminous BHs at 2; = 6 and 7. This 
suggests that the manner in which BHs populate halos is 
approximately redshift independent, save for the luminous 
BHs at redshifts 6 and 7. At these redshifts, we find the 
rarest (i.e. brightest) objects to be brighter than our 2 = 5 
fit would expect. In other words, at 2 = 6 and 7, the BHs 
in the most massive halos are more luminous than those in 
comparable halos at other redshifts, a result of unusually 
high accretion rates. To quantify this discrepancy, we per- 
formed HOD fits for 2 = 6, 7 (Table [2)| which are shown as 
dotted lines in Figure [T] We note that A^ remains approx- 
imately constant, suggesting that BHs of ~ IO^'^Lq tend 
to populate similar-size halos regardless of redshift (at least 
for 2^5). However, at 2 = 6 and 7 the smaller Bf, implies 
that BHs found in massive halos at 2 = 6 and 7 tend to be 
brighter than those in comparable halos at other redshifts, 
and the effect gets stronger for higher mass halos/higher 
luminosity BHs. 

This can be expl ained by consideri ng the growth history 
of the massive BHs. iDi Matteo et al.l (2011) find that the 
most massive BHs typically undergo a period of rapid (i.e. 
Eddington) growth in this general redshift range (2 ~ 6 — 8) 
as a result of high density streams of cold gas which not 
only help assemble the first massive halos, but appear to 



(1) 
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Figure 1. Left panel: Quasar lumino sity function for z = 11 to z = 5 (solid lines wi th Poisson error bars), along with the z = 5 and 6 
observational datapoints compiled by IWillott et alj ll2010al ) and iHopkins et alj 1120071) [See text for complete references to observational 
data]. We also show the predicted QLF from our HOD model using the z = 5 fit (dashed lines) and using redshift-specific fits (dotted 
lines). [See Table |2]for fitted parameters.] Right panel: Quasar magnitude function for z = 11 to z = 5, based on apparent bolometric 
magnitude. The regions probed by CANDELS and VIDEO surveys are shown as shaded regions (solid red - CANDELS DEEP; hatched 
blue - CANDELS WIDE; hatched green - VIDEO). The lower limits are based on the survey volume for z it 0.5. We also show the 
magnitude limits for JWST NIRCam (black dashed), LSST (pink dot-dashed). Dark Energy Survey (red long dashed), and VIKING 
(yellow triple-dot dash). 



penetrate to the halo centers. These cold flows facilitate ex- 
tended periods of Eddington growth of BHs in the most 
massive halos, which manifests itself in the QLF by flatten- 
ing the bright end. This continues until the energy released 
by the quasar heats the surrounding gas to such a point that 
it is blown away fr om the halo ce nter, and the BHs become 
self-regulated (see iDi Matteo et al.l I2OIII . for an investiga- 
tion into this behavior), at which point the bright-end will 
steepen once again. In this picture, the QLF should roughly 
follow the evolution of the halo mass function except for 
those BHs undergoing their Eddington growth phase, where 
the luminosity should be unusually high. We find exactly 
these results in Figure [ll where the z = 5 fit does a good 
job except in the bright-end at 2 = 6 — 7, where Eddington 
growth is common. 

In addition to the luminosity function shown in the 
left panel of Figure [l] we show the magnitude function in 
terms of apparent bolometric magnitude in the right panel 
which we use as a basis for predictions for several upcom- 
ing surveys. In particular, we show the regimes probed by 
the Cosmic Assembly Near-IR Deep Extragalactic Legacy 
Survey (CANDELS) WIDE and DEEP surveytQ (solid red 
and hatched blue, respectively). These regions are bounded 
by the magnitude limit of the F125W filter (converted to 
bolometric m agnitude using the spectral energy distribu- 
tion (SED) of lHopkins et al.ll2007l ). and the volume enclosed 
by the survey areas (m < 26.4 over 0.2 deg^ for WIDE 
and m < 27 A over 0.04 deg'^ for DEEP, assuming a red- 
shift range of ±0.5). We predict that both the WIDE and 
DEEP programs will find sources out to z ~ 7 — 8, with 
each probing a slightly different region of the luminosity 
function. We expect there to be ~ 119 ± 28 (~ 87 ± 28) 
quasars bright enough in the F125W band of the CAN- 
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DELS WIDE (DEEP) fields from z = 5-6, ~19±7 
(~ 18 ± 9) from 2 = 6 - 7, and ~ 1.7 ± 1.5 (~ 1.5 ± 1.5) 
from 2 = 7 — 8. This will drastically increase the number 
of high-redshift quasars, which will substantially improve 
the measurements of the faint end of the 2 = 6 QLF, and 
the different (though overlapping) ranges probed by the 
WIDE and DEEP programs will help constrain the faint- 
end slope. We have also provided a similar bounded vol- 
ume for the VISTA Deep Extragalactic Observations Survey 
(VIDEO), and magnitude limits for JWST's Near-Infrared 
Camera (NIRCam; black dashed line) (Gardner et al. 200q), 
Large Synoptic Survey Telescope (LSST; pink dot-dashed), 
Dark Energy Survey (red long dashed), and VISTA Kilo- 
Degree IR Galaxy Survey (VIKING; yellow triple-dot dash), 
each conv erted to bolometric magnitudes (using the SED of 
iHopkins et al. 2007 ) to provide the means for predicting the 
number of high-z quasars these surveys should find. Overall, 
these surveys should drastically improve the measurements 
of the faint-end QLF by improving the number of observed 
quasars, and pushing to luminosities an order of magnitude 
lower than current observations. 



3.2 Quasar Clustering 

In addition to the number density of quasars, we look at the 
clustering properties characterized by the correlation length 
(ro, defined as the scale at which the correlation function 
^(ro) = 1). In Figure [5] we show the predicted correlation 
length as a function of lower- luminosity cut for 2 = 1 to 6. 
These predicte d correlation lengths are obtained using a halo 
mas s functio n (^Rccd ct al. 2007) and the halo bias factor of 
Shet h et al.l (12001), weighted by the fraction of halos hosting 
sufficiently luminous quasars (Eqn. [ij using the 2 = 5 pa- 
rameters for 2^5 and the 2 = 6 fit at 2 = 6). We note that 
the halo bias facto r is not well constrained for the very high - 
mass end (see, e.g. [Tinker et al.ll2010l : IPillepich et al.ll2010l ). 



© 20?? RAS, MNRAS 000, ??-?? 



Early BHs in Cosmological Simulations 5 




Figure 2. Main panel: The predicted correlation length based on 
our HOD model as a function of luminosity (dashed lines). Curves 
of constant apparent i-band magnitude are shown in black, and 
the magnitude cut used for quasars by SDSS (19.1 for 2 < 3, 
20.2 for 2 > 3) is shown in red. Note that the z < 5 curves 
use the 2 = 5 HOD fit. Inset: Comparison of ro from the HOD- 
based prediction (dashed lines) and that obtained directly from 
our simulation (solid lines) for 2 = 5 (red) and 6 (orange). 
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Figure 3. Evolution in the predicted correlation length assum- 
ing a redshift independent HOD fit (using 2 = 5 parameters) for 
magnitude cuts of 19.1 (black, SDSS QSO cut for 2 < 3), 20.2 
(pink, SDSS QSO cut for 2 > 3), and 23.0 (blue). The dashed line 
from 2 = 5 to 6 shows ro if 2 = 6 HOD fit is used at 2 = 6. Filled 
circles: Observational data from lShen et al.l 112009 1 obtained with 
and without negative datapoins (dashed and solid error bars, re- 
spectively). 



and has not been investigated siRnifi cantly at the s e high red- 
shifts. Nonetheless, we find that the lSheth et ahl l|200ll ') bias 
factor does a good job of modehng our halo bias, though 
at 2: ~ 6 it underestimates by ~ 5% relative to our simu- 
lation (note that we adjustment the bias factor by 5% to 
account for this when calculating to^bh)- This technique is 
important as it allows us to extend our predictions to higher 
luminosities where we do not have sufficient BHs to obtain 
clustering statistics. 

To confirm the validity of this approach, we calculate 
the correlation length directly from the simulation using a 
maximum likelihood estimator (see, e.g. ICroft et al.lll997^ ■ 
and show the comparison in the inset plot (dashed lines - 
HOD model; solid lines - direct from simulation), finding 
good agreement in the luminosity range probed directly by 
our simulation. Thus we confirm that the large-scale clus- 
tering of luminous quasars can be well-modeled by the clus- 
tering of dark matter halos and a relatively simple BH HOD 
(Eqn. [1]). In the main plot we also show curves of constant 
apparent magnitude (solid black curves), and a curve for the 
magnitude cut used for quasar selection in SDSS {rrii < 20.2 
for z > 3, mi < 19.1 for z < 3, solid red curve). These curves 
show that using a magnitude cut rather than a luminosity 
cut will increase the observed redshift evolution. We also see 
that the evolution of the clustering strength with luminos- 
ity is relatively minor for Iow-Lbh (where BHs occupy the 
slowly-evolving end of the halo mass function) , but becomes 
more significant at high luminosities where BHs occupy the 
high-mass tail of the mass function. 

In addition, we note the effect of the change in the HOD 
model from z = 5 (red) to z = 6 (orange) . For low luminosi- 
ties (below ~ IO^^Lq) there is minimal effect, due to the 
similarity in A^. However, the smaller B^ at 2: = 6 (Ta- 
ble [2)1 means that higher luminosity BHs will typically be 
found in smaller halos, which have a correspondingly lower 



correlation length. Because the effect grows with luminosity, 
the difference manifests itself in the shallower slope of ro vs. 
Lbh, even resulting in a crossover at '^ 4 x 10^^ Lq. We note 
that this crossover is highly dependent on the HOD at high 
luminosities, and we are forced to extrapolate our HOD to 
reach these values. Thus this crossover is not a strong pre- 
diction (see further discussion below) but the general sup- 
pression is a clear trend in our simulation. In this way we see 
that ro could be effected by a change in the HOD, but the ef- 
fect should be fairly minor and will only occur at extremely 
high luminosities ( ^ 10^'^Lq). 

In Figure [3] we show the evolution in our predicted cor- 
relation length for SDSS-type quasars based on our z = 5 
HOD fit (for luminosity cuts corresponding to rrii < 19.1 
- Black; rm < 20.2 - Pink; rui < 23.0 - Blu e ), tog ether 
with the observed measurements of IShen et al.l ( 20091 ) (ob- 
tained with and without inclusion of negative datapoints, 
shown with dashed and solid error bars respectively). The 
curves assume constant occupation behavior (i.e. the z — 5 
HOD parameters are used at all redshifts). Based on ear- 
lier work, we would expect relatively minor evolut ion in the 
HOD between z — 5 and 2 = 3 (jChatteriee et al.ll201ll ). so 
the extrapolation should remain fairly accurate for 2 > 3. 
This extrapolation shows that our simulation is able to re- 
produce the large correlation lengths (ro ~ 17 /i~^ Mpc at 
2 = 4) from observational constraints. In addition, we match 
the rough evolution with redshift for 2 > 3, and we expect ro 
to continue to increase with redshift, reaching ~ 22h~^ Mpc 
by 2 = 6. Because we use a redshift-independent BH HOD 
model, this agreement with the observed evolution of ro with 
redshift suggests that high-z quasar clustering evolution can 
be completely explained by the evolution in the clustering of 
dark matter halos, i.e. without evolution in the mass of the 
typical host halo. Be low 2 = 3 we expect th e AGN HOD to 
evolve more quickly IjChatteriee et al.ll201ll ). so the extrap- 
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Table 2. Luminosity dependence of ^ and a (see Eqn.[TJ. 

fJ.,(7{L-£ 



(-^BH,cut 



Redshift 



Bu 



B„ 



5 27.05 

6 26.905 

7 27.07 



.0208 
.0191 
.0173 



.209 .128 
.239 .133 
.298 



olation should not be considered an accurate prediction at 
lower redshifts (the dotted curves of Fig.[3|, but we nonethe- 
less appear to remain broadly consistent with observations. 
In addition to the curves for the 2 = 5 HOD, we also 
show the effect of the change in occupation distribution from 
z = 5 to 6 (dashed lines). At z = 6 the shallower slope in 
the Z/BH — A/Host relation means that a given luminosity cut 
will correspond to a smaller typical host mass (and thereby 
produce a smaller ro) at 2 = 6 than 2 = 5. Because the 
difference in the IfOD is in the slope (B^j) of the -Lbh — Mhosi 
relation, the suppression gets stronger at higher luminosities 
and has essentially no effect on lower luminosities, which 
we show with the m; = 23.0 curve of Figure [S] We note 
that the magnitude of the suppression is highly dependent 
on the slope of ^ in Equation [T] (found in Table [2]) and is 
only significant at luminosities above those well-probed by 
our simulation. As such, we cannot accurately estimate the 
magnitude of the suppression (or the exact luminosity at 
which the 2 = 5 and 6 curves of Figure [2] will cross) , but 
we nonetheless expect a suppressed slope in the ro evolution 
from 2 = 5 to 6. Given the approximate magnitude of the 
suppression found in our simulation (on the order of ~ 5%) 
and the difficulties in measuring ro to high precision at such 
high redshifts, it is unlikely that such suppression will be 
detected in the near future, but it is nonetheless an effect 
which should be considered wherever possible. 



4 CONCLUSIONS 

In this paper we have investigated high redshift black holes 
found within a new large-scale hydrodynamic simulation, fo- 
cusing on the luminosity function and clustering behavior. 
We model the QLF for z ^ 5, probing black holes to lumi- 
nosities up to -^ 10^'* L©. We find reasonable agreement with 
observational data at z = 5 and 6, generally falling within 
the variation between surveys, confirming that our model is 
capable of producing Sloan-type quasars at very early times 
as well as matching their observed abundances. 

Using a HOD fit for the central AGN at z = 5 (which 
we have provided), we find that the evolution in the QLF 
is well described by the evolution in the halo mass func- 
tion (at least for the redshifts modeled by our simulation) 
without significant evolution in the occupation distribution, 
except for the high-luminosity end at z = 6 — 7 where we 
find a significantly flatter luminosity function. We postu- 
late this flattening of the luminosity function to be a re- 
sult of the largest black holes tending to undergo unusually 
rapid growth (and thereby producing unusually high lumi- 
nosities) during these redshifts, a conclusion supported both 
by the difference in the HOD fits at these redshifts, and by 
direct investigation into the lightcurves of these black holes 



(|Di Matteo et al.ll201 j ). At lower redshift (by z ~ 5), self- 
regulation suppresses the most massive BH growth, resulting 
in fainter bright-end BHs, and a corresponding steepening 
of the QLF. In particular, we note that this increase in lumi- 
nosity at z ~ 6 — 7 should make observations of such high-z, 
luminous (> IO^^Lq) BHs easier than would otherwise be 
expected. 

We used our luminosity function to provide estimates 
on the number density of detectable high redshift quasars 
for several upcoming surveys. In particular, we expect the 
CANDELS survey to find quasars up to z ~ 7 in both the 
WIDE and DEEP programs, with each probing a slightly 
different region of the QLF. At z '^ 6 these programs should 
each detect dozens of quasars across a wide range of lumi- 
nosities, drastically improving the observational constraints 
on the faint end of the high redshift QLF, particularly the 
faint-end slope. In addition, we have provided our estimated 
number density at each redshift for several additional up- 
coming surveys, thereby providing the approximate number 
density of quasars to be found as well as the survey areas 
necessary to reach the highest redshifts. 

We also investigate the clustering behavior of these 
high redshift quasars, finding luminosity dependent corre- 
lation lengths on the order of ~ 10 — 14 ft~^Mpc. Using 
our HOD model and the theoretical clustering of dark mat- 
ter halos, we show the correlation length as a function of 
BH luminosity, finding relatively weak luminosity depen- 
dence at low Lbh, but with increasing LBH-dependence at 
higher luminosities (where black holes are typically found 
in halos in the steep end of the halo mass function). We 
also compare our predicted ro to high-redshift observations 
and find excellent agreement, with our simulation predict- 
ing ro ~ 15 — 20 h~^ Mpc for Sloan- type quasars. We 
also roughly match the evolution of ro with redshift using 
a redshift-independent HOD, suggesting that high-redshift 
quasar clustering evolution is fully explained solely by the 
evolution in clustering of dark matter halos without a change 
in the typical host halo mass. 

Finally, we note the effect a change in the BH occupa- 
tion distribution from z = 5 to 6 has on ro, finding that it 
should suppress the clustering strength of high luminosity 
quasars at z = 6. Our limited sample size and evolution in 
the halo bias factor found in our simulation make it diffi- 
cult to quantify the magnitude of this suppression, but we 
nonetheless conclude that some suppression should occur 
(though likely below the sensitivity of upcoming surveys). 



ACKNOWLEDGMENTS 

This work was supported by the National Science Foun- 
dation, NSF Petapps, OCL0749212 and NSF AST- 
1009781. The simulations were carried out on Kraken 
at the National Institute for Computational Sciences 
(http://www.nics.tennessee.edu/). 



REFERENCES 

Barger A. J., et al., 2003a, AJ, 126, 632 

© 20?? RAS, MNRAS 000, ??-?? 



Early BHs in Cosmological Simulations 7 



Barger A. J., Cowie L. L., Capak P., Alexander D. M., 

Bauer F. E., Brandt W. N., Garmire G. P., Hornschemeier 

A. E., 2003b, ApJL, 584, L61 
Barger A. J., Cowie L. L., Mushotzky R. F., Yang Y., Wang 

W.-H., Steffen A. T., Capak P., 2005, AJ, 129, 578 
Begelman M. C, Volonteri M., Rees M. J., 2006, MNRAS, 

370, 289 
Bondi H., 1952, MNRAS, 112, 195 
Bondi H., Hoyle F., 1944, MNRAS, 104, 273 
Bonoli S., Marulli F., Springel V., White S. D. M., Bran- 

chini E., Moscardini L., 2009, MNRAS, 606 
Booth C. M., Schaye J., 2009, MNRAS, 398, 53 
Bower R. G., Benson A. J., Malbon R., Helly J. C., Frenk 

C. S., Baugh C. M., Cole S., Lacey C. G., 2006, MNRAS, 

370, 645 
Bromm V., Larson R. B., 2004, ARA&A, 42, 79 
Bromm V., Loeb A., 2003, ApJ, 596, 34 
Brown M. J. I., et al., 2006, ApJ, 638, 88 
Burkert A., Silk J., 2001, ApJL, 554, L151 
Chatterjee S., DeGraf C, Richardson J., Zheng Z., Nagai 

D., Di Matteo T., 2011, ArXiv e-prints:1104.3550 
Churazov E., Sazonov S., Sunyaev R., Forman W., Jones 

C, Bohringer H., 2005, MNRAS, 363, L91 
Ciotti L., Ostriker J. P., 2007, ApJ, 665, 1038 
Cirasuolo M., Magliocchetti M., Celotti A., 2005, MNRAS, 

357, 1267 
Colberg J. M., di Matteo T., 2008, MNRAS, 387, 1163 
Cristiani S., et al., 2004, ApJL, 600, L119 
Croft R. A. C, Dalton G. B., Efstathiou G., Sutherland 

W. J., Maddox S. J., 1997, MNRAS, 291, 305 
Croft R. A. C, Di Matteo T., Springel V., Hernquist L., 

2009, MNRAS, 400, 43 
Groom S. M., et al., 2005, MNRAS, 356, 415 
Groom S. M., Smith R. J., Boyle B. J., Shanks T., Miller 

L., Outram P. J., Loaring N. S., 2004, MNRAS, 349, 1397 
Croton D. J., 2009, MNRAS, 394, 1109 
Croton D. J., et al., 2006, MNRAS, 365, 11 
da Angela J., et al., 2008, MNRAS, 383, 565 
DeGraf C., Di Matteo T., Springel V., 2010, MNRAS, 402, 

1927 
Degraf C, Di Matteo T., Springel V., 2011, MNRAS, 413, 

1383 
Di Matteo T., Colberg J., Springel V., Hernquist L., Sijacki 

D., 2008, ApJ, 676, 33 
Di Matteo T., Khandai N., DeGraf C, Feng Y., Croft R., 

Lopez J., Springel V., 2011, ApJL submitted 
Di Matteo T., Springel V., Hernquist L., 2005, Nature, 433, 

604 
Ebrero J., et al., 2009, A&A, 493, 55 
Escala A., Larson R. B., Coppi P. S., Mardones D., 2004, 

ApJ, 607, 765 
Fan X., et al., 2004, AJ, 128, 515 
Fan X., et al., 2001a, AJ, 122, 2833 
Fan X., et al., 2006, AJ, 132, 117 
Fan X., et al., 2001b, AJ, 121, 54 
Fan X., et al., 2003, AJ, 125, 1649 
Ferrarese L., Merritt D., 2000, ApJL, 539, L9 
Fiore F., et al., 2003, A&A, 409, 79 
Gardner J. P., et al., 2006, SSR, 123, 485 
Gebhardt K., et al., 2000, ApJL, 539, L13 
Graham A. W., Driver S. P., 2007, ApJ, 655, 77 



Granato G. L., De Zotti G., Silva L., Bressan A., Danese 
L., 2004, ApJ, 600, 580 

Hopkins P. F., Hernquist L., Cox T. J., Di Matteo T., Mar- 
tini P., Robertson B., Springel V., 2005a, ApJ, 630, 705 

Hopkins P. F., Hernquist L., Cox T. J., Di Matteo T., 
Robertson B., Springel V., 2005b, ApJ, 630, 716 

Hopkins P. F., Hernquist L., Cox T. J., Di Matteo T., 
Robertson B., Springel V., 2005c, ApJ, 632, 81 

Hopkins P. F., Hernquist L., Cox T. J., Di Matteo T., 
Robertson B., Springel V., 2006a, ApJS, 163, 1 

Hopkins P. F., Hernquist L., Cox T. J., Robertson B., Di 
Matteo T., Springel V., 2006b, ApJ, 639, 700 

Hopkins P. F., Hernquist L., Martini P., Cox T. J., Robert- 
son B., Di Matteo T., Springel V., 2005d, ApJL, 625, L71 

Hopkins P. F., Richards G. T., Hernquist L., 2007, ApJ, 
654, 731 

Hoyle F., Lyttleton R. A., 1939, in Proceedings of the 
Cambridge Philosophical Society, vol. 35 of Proceedings 
of the Cambridge Philosophical Society, 405 

Jiang L., et al., 2009, AJ, 138, 305 

Johansson P. H., Naab T., Burkert A., 2008, Astronomische 
Nachrichtcn, 329, 956 

Kauffmann G., Haehnelt M., 2000, MNRAS, 311, 576 

Kawata D., Gibson B. K., 2005, MNRAS, 358, L16 

Kormendy J., Richstone D., 1995, ARA&A, 33, 581 

La Franca F., Andreani P., Cristiani S., 1998, ApJ, 497, 
529 

La Franca F., et al., 2005, ApJ, 635, 864 

La Franca F., et al., 2002, ApJ, 570, 100 

Lemson G., Virgo Consortium t., 2006, ArXiv Astrophysics 
e-prints: 0608019 

Lidz A., Hopkins P. F., Cox T. J., Hernquist L., Robertson 
B., 2006, ApJ, 641, 41 

Lopez J., Degraf C, DiMatteo T., Fu B., Fink E., Gibson 
G., 2011, in Statistical and Scientific Databases Manage- 
ment Conference (SSDBM), Portland, OR 

Magorrian J., et al., 1998, AJ, 115, 2285 

Makino J., Funato Y., 2004, ApJ, 602, 93 

Malbon R. K., Baugh C. M., Frenk C. S., Lacey C. G., 
2007, MNRAS, 382, 1394 

MaruUi F., Bonoli S., Branchini E., Gilli R., Moscardini L., 
Springel V., 2009, MNRAS, 396, 1404 

Marulli F., Bonoli S., Branchini E., Moscardini L., Springel 
v., 2008, MNRAS, 385, 1846 

Mayer L., Kazantzidis S., Madau P., Colpi M., Quinn T., 
Wadsley J., 2007, Science, 316, 1874 

Myers A. D., Brunner R. J., Nichol R. C, Richards G. T., 
Schneider D. P., Bahcall N. A., 2007, ApJ, 658, 85 

Pillepich A., Porciani C., Hahn O., 2010, MNRAS, 402, 191 

Porciani C, Magliocchetti M., Norberg P., 2004, MNRAS, 
355, 1010 

Reed D. S., Bower R., Frenk C. S., Jenkins A., Theuns T., 
2007, MNRAS, 374, 2 

Richards G. T., et al., 2005, MNRAS, 360, 839 

Richards G. T., et al., 2006, AJ, 131, 2766 

Ross N. P., et al., 2009, ApJ, 697, 1634 

Sazonov S. Y., Ostriker J. P., Sunyaev R. A., 2004, MN- 
RAS, 347, 144 

Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337 

Shen Y., et al., 2007, AJ, 133, 2222 

Shen Y., et al., 2009, ApJ, 697, 1656 

Sheth R. K., Mo H. J., Tormen G., 2001, MNRAS, 323, 1 



© 20?? RAS, MNRAS 000, ??-?? 



B Colin Degraf et al. 

Sijacki D., Springel V., di Matteo T., Hernquist L., 2007, 

MNRAS, 380, 877 
Sijacki D., Springel V., Haehnelt M. G., 2009, MNRAS, 

400, 100 
Silverman J. D., et al., 2005, ApJ, 624, 630 
Silverman J. D., et al., 2008, ApJ, 679, 118 
Springel V., 2005, MNRAS, 364, 1105 
Springel V., Di Matteo T., Hernquist L., 2005a, MNRAS, 

361, 776 
Springel V., Hernquist L., 2003, MNRAS, 339, 289 
Springel V., et al., 2005b, Nature, 435, 629 
Tinker J. L., Robertson B. E., Kravtsov A. V., Klypin A., 

Warren M. S., Yepes G., Gottlober S., 2010, ApJ, 724, 

878 
Tremaine S., et al., 2002, ApJ, 574, 740 
Ueda Y., Akiyama M., Ohta K., Miyaji T., 2003, ApJ, 598, 

886 
Volonteri M., Haardt F., Madau R, 2003, ApJ, 582, 559 
Volonteri M., Rees M. J., 2006, ApJ, 650, 669 
Willott C. J., et al., 2010a, AJ, 140, 546 
Willott C. J., et al., 2010b, AJ, 139, 906 
Wyithe J. S. B., Loeb A., 2003, ApJ, 595, 614 
Yencho B., Barger A. J., Trouille L., Winter L. M., 2009, 

ApJ, 698, 380 
Yoshida N., Omukai K., Hernquist L., Abel T., 2006, ApJ, 

652, 6 
Zheng Z., et al., 2005, ApJ, 633, 791 



© 20?? RAS, MNRAS 000, ??-?? 



